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ABSTRACT 

This  article  opts  for  the  use  of  direct  simulation 
Monte-Carlo  methods  to  solve  such  problems  as  those  associated 
with  winding  flows  around  flat  plates  of  finite  length.  The 
methods  in  guestion  are  ones  which  go  through  tracking  the 
movements  of  simulated  molecules  by  computer  in  order  to  realize 
numerical  value  simulations.  Collision  calculations  between 
simulated  molecules  are  determined  from  statistical  samplings. 

As  far  as  collision  models  are  concerned,  respective  selections 
were  made  for  the  use  of  hard  sphere  molecule  models  and  inverse 
power  law  models.  As  far  as  simulation  molecule  and  solid  wall 
effects  are  concerned,  option  was  made  for  the  use  of  models 
which  are  made  up  from  mixtures  on  the  basis  of  ratios  of  fully 
diffused  reflection  and  specular  or  mirror  surface  reflection. 

For  the  sake  of  the  reliability  of  empirical  test  methods, 
calculations  were  also  done  of  the  one  dimensional  flow  movements 
associated  with  such  questions  and  Rayleigh  problems  and  shock 
wave  structures.  Two  dimensional  calculations  opted  for  the  use 
of  synchronous  parallel  programs.  Numerical  value  results 
clearly  showed  that  direct  simulation  Monte-Carlo  methods  were 
capable  of  relatively  good  simulations  of  a  number  of  problems  in 
rarefied  gas  mechanics .  As  far  as  two  dimensional  calculations 
are  concerned,  the  machine  time  expended  and  the  necessary 
content  are  both  within  ranges  permitted  by  domestic  machines. 


SYMBOLS 


L*  Flat  Plate  Length 

L  L*/Xm  Nondimensional 

Flat  Plate  Length 

m*  Degenerate  Mass 


Vp  Most  Probable  Velocity  in 
Undisturbed  Gases 

Vr  Relative  Velocity 

Wm  Nondimensional  Collision 
Parameter  Cut  Off  Value 


*  Numbers  in  margins  indicate  foreign  pagination. 
Commas  in  numbers  indicate  decimals. 
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M  Mach  Number  or  the 
Oncoming  Flow 

n  Simulation  Molecule 
Number  Density 

Nc  Simulation  Molecule 
Number  in  Grid 

p  Pressure 

t  Time 

T  Temperature 

T0  Overall  Temperature 
of  Oncoming  Flow 

u  x  direction  speed 

D  y  direction  speed 


p  Density 

A  Average  Free  Travel 
v  Inverse  Power  Law  Index 


Subscripts 

1  Value  In  Front  of  Wave 

2  Value  Behind  Wave 

fr  Free  Molecular  Flow  Full 
Diffusion  Reflection  Value 

w  Wall  Surface  Value 

co  Oncoming  Flow  Value 


INTRODUCTION 

Normally,  it  is  necessary  to  opt  for  the  use  of  Boltzmann 
eguation  descriptions  of  the  phenomena  produced  by  rarefied 
gases.  This  eguation  is  a  nonlinear  integral  differential 
equation  of  distribution  functions  related  to  seven  independent 
variables.  Besides  extremely  particular  special  cases,  equation 
analyses  and  solutions  are  extremely  difficult.  Moreover,  as  far 
as  the  application  of  finite  differential  methods  of  numerical 
value  solutions  is  concerned,  it  is  necessary  to  expend 
extremely  large  amounts  of  machine  time  as  well  as  requiring 
extremely  large  internal  storage.  Because  of  this,  it  is  not 
economical.  /325 

Looking  from  the  point  of  view  of  basic  physical  substance 
or  quality,  the  microprocesses  produced  in  rarefied  gases  are  all 
random  processes.  At  the  end  of  the  1950' s,  people  began  to  pay 
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attention  to  the  similarity  between  this  basic  physical  nature 
and  Monte-Carlo  simulation  methods.  Applying  Monte-Carlo 
experimental  particle  methods  to  simulate  rarefied  gas  dynamics 
problems,  results  were  obtained.  The  first  were  Alder  and 
Wainwright  applying  this  method  to  study  the  process  of  gases 
tending  toward  equilibrium  states  in  closed  systems.  After  that, 
Bird  and  others  made  very  great  improvements  and  developments  on 
Alder  methods,  setting  up  direct  flow  field  simulation  methods 
and  successfully  carrying  out  winding  flow  calculations  for 
complex  rotationally  generated  bodies  [1]  as  well  as  techniques 
for  the  display  of  simple  physical  shape  winding  flows, 
relaxation  processes,  shock  wave  structures,  shock  wave 
reflection,  and  rarefied  gas  flow  fields. 

This  article  opts  for  the  use  of  direct  simulation  Monte- 
Carlo  methods  to  carry  out  numerical  value  solutions  in  finite 
length  flat  plate  winding  flow  problems.  For  the  sake  of  the 
accuracy  of  the  testing  methods,  calculations  were  also  done  of 
such  one  dimensional  questions  as  Rayleigh  problems  and  shock 
wave  structures.  In  conjunction  with  this,  comparisons  were  made 
with  reference  results  and  theory. 


I.  SIMULATION  TECHNIQUES 

In  flow  fields,  a  control  volume  ADFC  was  selected  (see 
Fig.l).  It  is  assumed  that,  within  a  fixed  time  interval  Atm, 
there  are  AN  individual  molecules  which  cross  over  the 
boundary  surface  AD  and  enter  into  the  control  volume.  There  are 
AN'  individual  molecules  which  cross  over  the  boundary  surface 

CF  and  leave  the  control  volume.  In  the  control  volume  ,  N(t) 

A  r - nC 

I  I 

Hl_  i  ! 

i  i 

Fig.l  Control  Volume  1  ! 


individual  molecules  move  within  it  in  conjunction  with  which 
collisions  are  given  rise  to.  Monte-Carlo  direct  simulation 
techniques  are  to  go  through  and  carry  out  tracking  calculations 
for  these  molecules  and  collisions  between  particles.  Subsequent 
to  this,  statistical  averaging  obtains  the  relevent  macro 
quantities .  Due  to  machine  internal  storage  and  computer  time 
limitations,  it  is  not  possible,  nor  is  there  any  need,  to 
directly  track  all  actual  molecules  and  carry  out  calculations  in 
conjunction  with  that.  During  actual  simulations,  use  is  made  of 
simulation  molecules.  Each  simulation  molecule  represents  many 
individual  actual  molecules.  Due  to  this  being  the  case,  the 
number  of  molecules  in  the  control  volume  is  vastly  reduced. 
However,  there  is  a  normal  requirement  for  a  number  of  simulated 
molecules  in  an  individual  grid  which  is  not  less  than  20.  In 
another  regard,  within  the  time  interval  Atm  between  going 
from  t  to  Atm,  simulation  molecules  aperiodically  are  all 
capable  of  meeting  up  with  other  simulation  molecules  or  in 
collisions  with  solid  walls  and  altering  the  mode  of  their 
movements.  In  order  to  realize  simulation  calculations,  one 
assumes  that  simulation  molecules,  in  this  interval  of  time,  do 
not  give  rise  to  any  type  of  collision.  Moreover,  one  takes  the 
collisions  which  should  take  place  in  association  with  molecular 
systems  in  the  time  interval  Atm  and  delays  all  of  them  to 
completion  at  the  time  instant  t  +  Atm.  In  this  way,  when 
selecting  the  time  interval  Atm,  it  is  necessary  to  satisfy 
the  two  limiting  conditions  of  convection  and  collision  below 


U  A tJLx«  1  ,  V  A tJL,«  1 

Moreover,  one  must  also  satisfy 

Af.«*« 

In  these,  U  and  V  are,  respectively,  the  average  speed  components 
associated  with  the  x  and  y  directions.  Lx  and  Ly  are, 
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respectively,  the  characteristic  lengths  associated  with  the  x 
and  y  directions .  r'm  is  the  average  free  time  between  two 
instances  of  collision  between  simulation  molecules. 

1.  GRID  DIVISIONS  AND  INTIAL  FIELD  DEFINITION 

Statistical  averaging  and  collision  calculations  all  are 
carried  out  within  each  grid.  Because  this  is  the  case,  grid 
dimensions  cannot  be  too  large.  Otherwise,  the  macro  quantities 
obtained  from  statistical  averaging  would  be  unable  to  reflect 
macro  physical  configurations.  In  conjunction  with  that,  it 
wpuld  cause  a  loss  of  reality  for  speeds  from  distrubances  to 
propagation.  In  general,  the  required  grid  dimensions  Ax  and 
Ay  are  both  not  greater  than  the  average  free  travel  L 
associated  with  molecules  of  the  oncoming  flow.  This  article,  in 
calculations,  has  grid  dimensions  which,  by  contrast,  are  based 
on  determinations  associated  with  computer  internal  storage. 

Ax  and  Ay  are  selected  as  0.5  or  A,,,.  In  flat  plate  winding 
flow  calculations,  the  y  direction  overall  length  is  selected  as 
8-40  A*,.  The  x  direction  overall  length  is  selected  as 

50  -  100  A*,. 

/326 

At  a  zero  time  instant,  in  each  grid,  the  same  number  of 
simulation  molecules  are  positioned.  In  each  individual  grid, 
the  actual  location  of  simulation  molecules  is  determined  on  the 
basis  of  uniform  distribution  random  sampling.  The  speeds  which 
they  possess  are  determined  on  the  basis  of  equilibrium  state 
speed  distribution  function  sampling.  In  the  same  way,  when 
implementing  calculations,  the  number  of  simulation  molecules  in 
each  grid  must  be  determined  on  the  basis  of  computer  internal 
storage  and  the  level  of  complexity  of  the  problem.  In  carrying 
out  shock  wave  structure  problem  calculations,  at  zero  time 
instants,  shock  waves  are  positioned  at  the  exact  center  of  the 
control  volume.  The  thickness  is  zero. 
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2.  MOVEMENTS  OF  SIMULATION  MOLECULES  AND  COLLISIONS 
BETWEEN  SIMULATION  MOLECULES 

Within  time  interval  Atm,  the  movements  of  each 
simulation  molecule  maintain  uniform  speed  in  a  straight  line 

*  (/  +  £/„>=*  </)+«</) 
y  (t  +  Atj=  y  (/)  +  *;(/) A/. 

In  the  instant  t  +  Atm,  collision  calculations  are 
carried  out  within  each  grid.  Within  an  individual  grid,  random 
sampling  takes  out  two  simulation  molecules.  On  the  basis  of 
methods  for  acceptance  or  rejection,  it  is  determined  whether  or 
not  they  give  rise  to  collisions.  If  collision  occurs,  then,  on 
the  basis  of  the  conservation  of  mass,  the  conservation  of 
momentum,  the  conservation  of  energy,  as  well  as  the  collision 
deviation  angles  given  out  by  sampling,  according  to  two 
dimensional  collision  theory,  calculations  give  the  speeds 
associated  with  the  two  simulation  molecules  after  collision. 
Subsequent  to  this,  calculations  give  "collision  time"  for  the 
instance  in  question. 

At  =  -~{nWl{Klm*ynr-l)  • 


Here,  k  is  the  molecular  interaction  law  constant. 

For  each  grid,  one  time  measurement  device  was  set  up.  As 
far  as  each  instance  of  the  occurrence  of  collisions  is 
concerned,  one  then,  on  the  corresponding  time  measurement 
device,  adds  on  the  instance  of  "collision  time"  in  question. 
When  the  values  of  time  measurement  devices  exceed  t  +  Atm, 
within  the  time  interval  Atm,  all  calculations  are  completed 
for  the  collisions  which  should  occur. 
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In  the  calculations  of  this  article,  in  all  cases,  option  is 
made  for  the  use  of  hard  sphere  collision  models  and  inverse 
power  law  collision  models. 

3 .  UPPER  FLOW  BOUNDARY  TREATMENT 

On  the  basis  of  numerical  flux  density,  calculations  are 
carried  out  of  the  number  of  molecules  to  enter  the  control 
volume  within  time  interval  Atm.  Generally  speaking,  when 
actual  molecular  numbers  are  turned  into  simulation  molecular 
numbers,  the  values  obtained  from  calculations  are  not  whole 
numbers.  In  the  process  of  taking  calculation  values  and  turning 
them  into  integers  or  whole  number  values,  the  rounding  off  of 
fractional  portions  is  determined  by  acceptance  or  rejection 
methods.  Each  individual  simulation  molecule's  entry  time  into 
the  control  volume  and  its  location  for  penetrating  the  boundary 
are  determined  from  sampling.  The  speed  of  each  simulation 
molecule  is  determined  on  the  basis  of  sampling  of  average  speeds 
for  unidirectional  equilibrium  state  distributions  associated 
with  . 


4 .  LOWER  FLOW  BOUNDARY  CONDITION  TREATMENT 

Behind  the  lower  flow  boundary,  as  far  as  the  placement  of 
the  image  of  a  control  volume  simulation  molecule  system  is 
concerned,  all  virtual  molecule  speeds  and  speeds  of  the  original 
simulation  molecules  are  the  same.  If,  in  the  time  interval  Atm, 
the  original  molecules  have  escaped  from  the  lower  flow 
boundaries,  then,  the  serial  numbers  are  canceled.  If,  in  the 
time  interval  Atm,  there  is  a  virtual  or  image  molecule 
which  enters  the  control  volume,  then,  one  takes  that  virtual  or 
image  molecule  and  sees  it  as  a  real  simulation  molecule  and 
gives  it  a  registration. 
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5  . 


OTHER  BOUNDARY  TREATMENTS 


Exterior  boundary  conditions  are  lens/mirror  surface  or 
specular  reflection  conditions.  As  far  as  the  number  of 
simulation  electrons  which  escape  from  boundaries  is  concerned, 
adjustments  are  carried  out  on  the  basis  of  the  law  of  the 
conservation  of  the  total  number  of  molecules  within  the  entire 
control  volume. 

Solid  wall  conditions  are  made  up  on  the  basis  of  a  mixture 
according  to  the  ratio  between  the  whole  diffusion  reflection  and 
l^ns/mirror  surface  or  specular  reflection. 


II.  ANALYSIS  OF  RESULTS 


/327 

For  Rayleigh  problem  calculation  results  see  Fig.’s  2-4. 
These  are,  respectively,  gas  speed  along  the  y  direction, 
density,  and  temperature  distributions.  At  the  instant  t  =  0, 
flat  plates  use  the  speed  Cr;,  to  suddenly  start  up.  At  the 
instant  t  =  4.43125  (-^■/  VQ),  viscous  diffusion  has  already 
reached  a  location  approximately  6  average  free  travels  from  the 
surface  of  the  plate.  In  the  vicinity  of  the  flat  plate,  gas 
flow  speeds  are  small  compared  to  flat  plate  speeds.  Gas  flow 
temperatures  are  high  as  compared  to  flat  plate  temperatures. 

This  clearly  shows  that  gas  flows  on  plate  surfaces  show  the 
appearance  of  speed  slipping  and  temperature  jumps.  In  the 
Fig.'s  it  is  also  shown  that  there  is  complete  agreement  with  the 
results  from  Reference  [1]. 
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Fig. 2  The  Distribution  for 
Tangential  Velocity, 
Density,  and  Temperature 

Key:  (1)  Reference 

[1]  Values  (2)  Calculated 

Values  from  This  Article 


Fig. 3  The  Density  Distribution 

Key:  (1)  Refernce  [2]  Values  (2) 
Calculated  Values  from  This  Article 
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Fig. 4  The  Temperature  Distribution 

(1)  Flow  Direction  Temperature 

(2)  Transverse  Direction  Temperature 

( 3 )  Sum  Temperature 


Fig. 5  The  Density  Distribution 
Key:  (1)  Hard  Sphere  Model 


Calculation  results  for  shock  wave  structure  problems  are 
given  in  Fig. 5.  In  the  Fig.'s,  there  are,  respectively,  given 
gas  density  distribution  along  the  direction  of  flow  as  well  as 
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flow  direction  and  transverse  direction  temperature 
distributions .  The  calculated  curves  and  the  results  from 
Reference  [2]  agree  with  each  other.  The  asymmetrical  natures  of 
the  calculated  curves  are  given  rise  to  due  to  the  fact  of 
relaxation  processes  during  momentum  and  energy  transfers.  With 
regard  to  calculation  results  for  different  inverse  power  law 
indices  v,  they  clearly  verified  the  formula  for  approximate 
calculations  with  regard  to  the  influence  of  molecular  models  on 
shock  wave  thicknesses  given  in  Reference  [2] 

LIJLaMV-1' 

Here,  L  is  the  shock  wave  thickness.  Xa  is  the  average  free 
travel  in  undisturbed  gases.  Ms  (unclear)  is  the  shock  wave  Mach 
number . 

Winding  flow  calculation  results  for  flat  plates  of  finite 
length  are  shown  in  Fig.'s  6-8.  Fig. 6  gives  wall  surface 
pressure  distributions.  In  conjunction  with  this,  comparisons 
are  carried  out  with  experimental  results  from  Becker  [3].  This 
article  opts  for  the  use  of  simple  single  atom  and  molecule  gas 
models  to  carry  out  numerical  value  simulations.  There  is  mutual 
consistency  with  Becker's  experiments  with  the  inert  gas  Ar  as 
medium.  However,  in  order  to  show  the  influences  of  internally 
free  molecules,  in  the  same  Fig.,  it  gives  experimental  results 
taking  N2  as  the  medium.  From  the  Fig.,  it  is  possible  to  see 
that  calculations  for  15%  specular  reflection  agree  relatively 
well  with  experimental  results.  It  points  out  our  properly 
selecting  diffusion  reflection  and  specular  reflection  /328 

percentages.  There  is  a  great  correlation  with  the  degree  of 
accuracy  in  the  case  of  simulation  analogs.  On  the  rear  edges  of 
flat  plates,  due  to  this  article's  opting  for  the  selection  of 
flat  plate  finite  length  conditions,  in  making  a  comparision  with 
the  results  of  experiments,  there  is  a  relatively  low  pressure 
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Fig. 6  The  Surface  Pressure  Distribution 

Key:  1)  Full  Diffusion  Reflection  (2)  15%  Specular  Reflection 


value.  This  point  is  also  capable,  from  physics,  of  being 
explained.  Due  to  the  existence  of  flat  plates,  molecules  close 
to  the  surface  of  the  object  gather  together  causing  flat  plate 
surface  pressures  to  rise.  This  leads  to  experimental  values 
being  higher  than  calculated  values.  Reference  [3]  gives  the 
experimental  conditions  to  be,  for  Ar  gas,  =  12.66, 

=0.08  Tq .  This  is  consistent  with  this  article's  calculated 
conditions  of  M.  =  12.9,  Tw  =  0.08  T0.  For  N2  gas,  Mro  =  11.18, 

Tw  =  0.1  T0.  There  is  some  discrepency  here  with  the  calculated 
conditions  of  this  article. 


Fig. 7  The  Surface  Pressure  Distribution 
Key:  (1)  Full  Diffusion  Reflection 
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Fig. 8  The  Surface  Pressure  Distribution 
Key:  (1)  Full  Diffusion  Reflection 


Fig. 7  gives  a  comparison  of  calculation  results  under 
conditions  of  different  plate  lengths  L.  From  the  calculation 
results,  one  can  see  the  influence  of  plate  length  on  wall 
surface  pressures.  This  is  primarily  concentrated  on  the  rear 
part  of  the  plate.  On  the  most  forward  portion  of  the  plate, 
wall  surface  pressures  receive  no  influence  from  plate  length  at 
all. 


What  is  worth  paying  attention  to  is  the  fact  that,  besides 
trailing  edges,  the  nondimensional  pressure  values  for  the  whole 
plate  are  all  greater  than  1.  This  clearly  shows  that  collisions 
between  molecules  as  well  as  between  molecules  and  wall  surfaces 
cause  wall  surface  pressures  to  rise. 
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